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Abstract 

This paper presents the very large eddy simulations (VLES) of a Jet-A spray reacting flow in a single 
element lean direct injection (LDI) injector by using the National Combustion Code (NCC) with and 
without invoking the Eulerian scalar DWFDF method, in which DWFDF is defined as the density 
weighted time filtered fine grained probability density function. The flow field is calculated by using the 
time filtered compressible Navier-Stokes equations (TFNS) with nonlinear subscale turbulence models, 
and when the Eulerian scalar DWFDF method is invoked, the energy and species mass fractions are 
calculated by solving the equation of DWFDF. A nonlinear subscale model for closing the convection 
term of the Eulerian scalar DWFDF equation is used and will be briefly described in this paper. Detailed 
comparisons between the results and available experimental data are carried out. Some positive findings 
of invoking the Eulerian scalar DWFDF method in both improving the simulation quality and maintaining 
economic computing cost are observed. 

1.0 Introduction 

Many engineering applications of large eddy simulations (LES) for internal combustor flows need to 
accurately account for the turbulence-chemistry interaction to facilitate a higher fidelity analysis of the 
design. In the conventional computational fluid dynamics (CFD) simulation methods the turbulence- 
chemistry interaction appears in an unclosed term that must be modeled, and its various empirical models 
have been proposed. The simplest one is the so called “laminar chemistry” which simply ignores the 
detailed effects of turbulence on the chemical reactions. However, in the DWFDF method, the chemical 
reaction source term in the DWFDF equation is in a closed form, therefore it can be directly calculated 
without any modeling, which is the unique feature of the probability density function (PDF) like 
approaches (Refs. 1, 2, and 3). 

In this paper, we present the preliminary VLES results of a complex Jet-A (C12H23) spray reacting 
flow in a single element LDI injector. The simulation methods include both the conventional CFD 
method (TFNS) and an Eulerian scalar DWFDF method (Ref. 4), in which the velocity field is determined 
by the continuity and momentum equations of time filtered compressible Navier-Stokes equations. All 
simulations are done with the NCC code (Ref. 5) using the same numerical parameter setting and the 
same computational domain and grid resolution. In the pure TFNS simulation, all the subscale turbulent 
flux models, i.e., stresses, heat and species fluxes, are nonlinear models (Refs. 6 and 7). The basic 
equations and models are described in Section 2.0. In the Eulerian scalar DWFDF method, the convection 
term of scalar DWFDF contains a conditional 'mean' that needs to be modeled. Here a nonlinear model is 
introduced, which provides a diffusion process in the sample variable space (Refs. 3 and 8). A brief 
description about the scalar DWFDF equation and the model is given in Section 3.0. 
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At the present stage of simulations, we use a relatively simple spray injection model (i.e., prescribed 
droplet size and distribution) and a global five-species one-step kinetics for combustion chemistry 
(Ref 9), so that we could concentrate more on the evaluation of different simulation approaches. 

The main results of simulations, including both with and without invoking an Eulerian scalar 
DWFDF method, are described in Section 4.0. They are compared side by side to examine how much 
effects are produced by invoking the Eulerian scalar DWFDF method. 

One of the main objectives of this study is to search for a consistent and stable VFES simulation tool 
for multiphase combustion flows. We require that this tool is of a rational “physics” based model and also 
able to produce the physically reasonable, numerically stable solution that can sustain over a very long 
(infinite) time period. And also, this tool should be relatively economic. It is encouraging to observe from 
the present simulations that the adopted Eulerian scalar DWFDF method does show the potential for 
improving the simulation quality and remaining an economic computing cost. 

2.0 Basic Equations for TFNS Simulations 

2.1 Time Filtered Turbulent Variables ((>( jc, r) ,<|> (x,t) 

In the case of compressible turbulent reacting flow, we often deal with two types of time filtering: one 
with the density weighting, the other without the density weighting. The filtered turbulent variable 
without the density weighting is denoted by ([)(-M) and is defined as 


— P+OO 

4>(jc, t ) = (|)(jc, t') G(t -t r )dt r 

J— 00 


( 1 ) 


where § is the unfiltered turbulent variable, e.g., velocity components U i9 density p, pressure P, species 

N 

mass fraction O t and internal energy e = . The integration is over the entire time domain 

m = 1 

- oo < f < + oo. G(t - f) is the time filter with a constant filter width A T and satisfies the following 
condition and asymptotic property: 


p+oo 

f G(t-t')dt' = 1 

J —oo 

( 2 ) 

r+oo 

§(x 9 t)G(t-t)dt =§(x 9 t) 9 asA r ^0 

J —00 

( 3 ) 


The density-weighted filtered turbulent variable is denoted by § (x,t) and is defined as 

(4) 

P 

These filtered variables (K*,0 , § (x,t) mainly represent the large scale turbulence, and they are still 
random but contain relatively low frequency part of the turbulent motion when comparing with the 
unfiltered turbulent variable ct>(jc,^) . 

2.2 Basic Equation 

Applying the time filtering, (1) and (4), to the compressible Navier-Stokes equations, we obtain: 


dp , dp Uj Q 

dt dxj 


( 5 ) 
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where 


dpUj | dpUjUj dP_ | 8 


dt dxj 


Ox.- dXj 

1 j 


2pv(i 9 --8,S a ) 


ep i+ apc if= _^ + — + 

dt 9x- 9x ; . kk 


r 






V 


+e 


y 


9pO m 0p[/-O m 9 


dt 


- + 


9x- 9x- 


p(m) ^O m 

& i J 


+ W m m = l,2,-*-,M 


p.p R ^ m T _pR^ m e 9 or p = E^°M + i 




( 6 ) 

(7) 

( 8 ) 


(9) 

( 10 ) 


In the above equations, k, v and r (m) are the molecular heat conductivity, kinematic viscosity and the m-th 
species diffusivity. It is commonly assumed that r (m) is same for all species O ffl , and w m is the molecular 
weight. The h m9 T are the enthalpy of species and the temperature, Q is the radiation rate, W m = pS m is the 
chemical production rate of the m-th species, O m+1 represents the internal energy e 9 R is the universal gas 
constant. Note that all the fuel-spray related terms are not shown in above equations for the simplicity, 
please see Reference 10 for the details. These equations are general; however, unlike the constant density 
flows, further approximations for the terms on their right hand side are required in order to complete the 
density weighted time filtering process. One of such approximations leads to 


2pv 



3 j 


11 


dUi 

dxj 


+ 


dUj_ 

9x- 


2- df/9 

— Off 

3 J dx k 


( 11 ) 


In which, we have basically neglected the variations of pand pduring the filtering process, the value of p 
will be considered as a known function of P 9 T, • • • . Other types of approximations are also possible, for 
example, 


2pv 


£, 8 a Si 


ij°kk 


"dp Ut , dpUj 2 § 

^ dxj 9x- 3 ij 


dp ud 

dx k y 


( 12 ) 


Similarly, 


% 





m - 1 


90 

^ m 

9x- 


k 9pT 

p dx < 


Z r< "p' , « 


m - 1 


90 

^ m 

9x- 


(13) 


pp(m) *^O m ~ p(m) ^P ^ m 

dx t dx i 


(14) 


Where v, k, c and T (w) are considered as the known functions of P 9 J 7 , • • • . Furthermore, we introduce the 
dissipation rate of the subscale turbulent kinetic energy, 
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Therefore, the time filtered compressible Navier-Stokes equations can be written as 


dp | dp Uj _ Q 
dt dxj 

dp Uj dp Up) dP d 
dt dxj 8x t dxj 


' dpU± | dp Uj 2 dpUk/ 
dxj dx i 3 ij dx k 


Ct OX: OX; 


dp 0m t dpt/,.<D m d 

dt dX; dX; 


p(m) ^pOm 


8X; 


+ p5* w 


m = 1,2,-**,M 


'* J 


pR ^ O m O M+1 


P = p^V = P*|^ ? or P = 


m=l W m C v m= 1 


C v m =l 


C U =- 


k dpT 
p dx : 


M 


-Zp rl ”’*. 


i m = 1 


so 

jji 

3X; 


( 16 ) 

(17) 

(18) 

(19) 

( 20 ) 
(21) 


These equations are considered quite general, because i) they are exact if the flow becomes a constant 
density incompressible flow, ii) all the approximations made in Equations (17), (18) and (19) are related 
only to the molecular diffusion terms which are less important and even become negligibly small, because 
they are of 0(1 /Re) comparing with the convection terms on the left hand side for turbulent flows at high 
Reynolds numbers (see Refs. 1 1 and 1). In addition, Equation (16) to (19) are often used together with the 
further approximations for (20) and (21): 


M 


p= p*Z 


O T 


pRT\ pRT _l 

v M ) M 


= V^l _L = V 

M 


0/7 


M — 


_ k dpT 

q i=-~— 

P dx - 


( 22 ) 

(23) 


The momentum flux pC/ z -C/ • , the energy flux pU t e and the species flux pU i O m are considered to be 

critically important and should be carefully modeled. Many models in the literature, from the simplest 
Smagorinsky model (Ref. 12) to the complex two-equation nonlinear models (Refs. 13 and 14), have been 
suggested in terms of subscale turbulent stresses and subscale turbulent scalar fluxes: 

t- = p (UjUj - ufij) and ©• = p(C/-0 - f?-0) , where 0 represents the scalar quantities: e and O m . 


2.3 Nonlinear Models for Subscale Turbulent Fluxes 

A general constitutive relationship between subscale turbulent stresses Xy and the strain rate of filtered 
turbulent flow Sy , Cl- suggests (Ref. 15) 
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( 24 ) 


%-{V*=-2/C m Py(^- 6 A/ 3 ) 

i 4 

+ 24 5 / p -^- [^5^- - Sf k Q k j + Ci. ik S km Q. mJ - Cl kl S lm Cl mk 8 jj + II s (Sy - 5 yS kk /3)J , 

where, 5, =(c/ fJ + C/,y)/2 , Cly=(Vij -U u )/2 , II s =(S kk ~S mm -~S kl ~S lk )/ 2 .Themodel 

coefficients C^, ^4 3 and ^4 5 are constrained by the realizability condition and the rapid distortion theory 
limit. They are formulated as (see Ref 14): 


C„ 


4.0 + A S -U 
8 


1.0-A 2 C 2 \ -S' 

b pi 


A, 


1.6Qp 


_r 


o.5 + i.5^-n*s* 

8 2 


J, 


5 —k 4 7S*S* + Q*Q* ’ 

P A — 4 — 


(25) 


in which, 


A s = V6coscp, cp = ^-arccos( 




TT/* _ °ij°jk^ki 

(sy ’ 

(26) 


(27) 


The coefficient /is a function of the resolution control parameter (RCP) that is defined as a ratio of the 
time filter width A T to a global integral time scale of the flow T: RCP = A T IT and 


/ 


f A 0 

1 « 2 1 

f A 0 


f Ar l 

T J 


l T ) 


l T J 


(28) 


As it is discussed in Reference 7 that RCP represents the percentage of unresolved turbulent kinetic 
energy. Therefore, the value of RCP and the coefficient / are always between 0 and 1 . Similarly, the 
nonlinear model for subscale scalar fluxes is formulated as (Ref. 16) 


0. - 0 r 


dp0 

dX: 


& T ( C Ay + c 2^ij) 


dp 9 


(29) 


Where denotes the subscale turbulent diffusivity for the corresponding scalar quantity 0. It is often 
approximated by S r = v/Pr 9 , and Pr e represents the turbulent Prandtl number or Schmidt number 
depending on whether the scalar quantity 0 is the energy e or the species . The subscale turbulent 

eddy viscosity is defined as v T =f • • k 2 /c. The coefficients, c x and c 2 , are yet to be calibrated. In the 

current simulations they are set to be c x = c 2 = -0.24. 

The subscale turbulent kinetic energy and its dissipation rate k , 8 will be determined from the 
following model equations: 


— p k + 
dt 


P 11; li = 

dx t 1 


d 

dx ■ 




V? 


-ps 


(30) 
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( 31 ) 


d _ d _ „ d 

— p 8 H p W ■ : 8 = 

8t dx i dx i 




‘ Qi T i7^i7 “ ~~ Q2 " 
81 J k 82 k 


where C e i and C e2 are the model coefficients. We have adopted the commonly used values of C e i = 1.45 and 
C e2 = 1 .92 in the present work while keeping in mind that here \a T is a subscale turbulent eddy viscosity. 


3.0 Basic Equation for Scalar DWFDF 

In this Section, we will use the fine grained probability density function (FG-PDF) to define the 
density weighted time filtered fine grained probability density function (DWFDF), then explore the 
relationship between scalar DWFDF and time filtered turbulent scalar variables. This will provide the 
basis for establishing the transport equation for the scalar DWFDF. 

3.1 Definition of Scalar DWFDF 

3.1.1 Fine Grained Probability Density Function (FG-PDF) for Scalars (v|/ ; x, t ) 

According to Pope’s definition (Ref. 1), the joint FG-PDF for turbulent velocity and scalars (i.e., 
compositions or species mass fractions, internal energy) can be written as 

= 8(U(x,t) - V)5(0(x,t) - v|/) 

3 M + 1 ( 32 ) 

i - 1 m = 1 


Its marginal FG-PDF for scalars is 


M + 1 

II (v ; x,t) = 8(©(*,0 - \\i) = F[ 8 ( a> «( JC »0 _ ) 

m = 1 


(33) 


where 5 denotes the delta function, U(x,t) is the turbulent (random) velocity vector (Uu U 2 , C/ 3 ), 0(x,f) is 
the turbulent (random) scalar array (O b 0 2? - • •, O m , O m +i), for example, M species mass fractions and one 
internal energy O m +i = e\ the xj denote the physical space variable (x u x 2 , x 3 ) and the time t , V= (V u 
V 2 ,V 3 ) and \| / = (\|/i, \|/ 2 , - - - , \| / M , \|/m+i) are the sample space variables for U(x,t) and O(x,0 , respectively. 

3.1.2 Scalar DWFDF F 0 (\|/; x,t) 

We define the following density weighted time filtered fine-grained probability density function for 
scalars as 


P+OO 

F 0 (y;x,t)= p(x,t') fl (v| i;x,t')G(t-t')dt' 

J —00 

= J + ”p(x,/')5(0(x,r)-v|/)G(/-/')A' 


( 34 ) 


Obviously, the scalar DWFDF, (i| /; x,f), is still a random quantity. And it satisfies the following 

“normalization” property: 
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| Fq (i|/ ; x, t) d\\i = 1 1 p(x,t')5(® -y)G(t -t')dt' d\\i = J p[x,t')G(t -t')dt' = p (35) 

-00 -00 
poo 

Note that here means that the integration is over the entire domain of the sample space \| /. 

J —00 

3.1.3 Relationship Between Scalar DWFDF and Large Scale Turbulent Scalar Variables 

With the definition of scalar DWFDF described in Equation (34), we can exactly deduce the density 
weighted time filtered scalar turbulent variables that are defined in Equations (1) and (4). For example, 

| v|/F 0 (v|/; x,t)d\\f = | | i|/p(x^')S(®(x^')-i|/)G(^-E)<iF<ii|/ 

= J 0(x,t , )p(x,t , )G(t-t')dt' (36) 

= p® = p(x,?)®(x,0 

Equation (36) indicates that the left hand side is an operation (®) that defines the density weighed time 
filtered turbulent variable p® : 


(®) = j \\i ; x,t)dy = p(x,t)Q> (x,t) (37) 

For a function W(Q(x,t)), it is easy to verify that 

/ \ /• +00 

(W (®)) = } %)^(i)» ; x,t)d^ = p{x,t)W{H>{x,t)) (38) 

Furthermore, we may consider the derivatives V® as a new random quantity and legitimately write 

(V®) = pV® (39) 


However, because of the variable density, the “operation” ( ) does not have the differential commute 
property, i.e., 

(V®)*V(®) (40) 

because p V® A V^p® j . 

We can also write the density weighted time filtered (G/®,) for the joint velocity and scalar variables as 
{UjO,) = j_ + V. F Ut o(V,\ |/; x,t)dVdy = p(j«;*)c7^(x,0 (41) 

where, (V 9 \|/; x,t) is the joint DWFDF defined as 


which also satisfies the “normalization” property: 

00 00 

J J Fu,<t>(V>W’ x,t)dVdy = j p(x,t')G(t -t')dt' = p (43) 

-00 —00 
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3.1.4 Joint DWFDF, F v & (F, \|/; x,t) and Its Conditional DWFDF $ (F|\|/; x,t) 

From Equation(42), we may follow Reference 1 to define a “conditional” DWFDF on the condition 

<D = \|/ as 


(^| 


Fjjd, ( V ,\|/; x,t) 
F 0 (v|/; x,t) 


(44) 


and the “conditional filtering” (or conditional ‘mean’) as 

(U(x, 0|v|/} = j_ + ” V F^ (F| v| ,;x,t)dV = J_ + ” F F u><Sf (V, y;x,t) dV 

, X,t ) 00 

= — J +C °F| +00 p(x,t')5(f/(x,t') -F)s(d>(jc,t') - X| i)G(t-t')dt'dV (45) 

= — J +CO p(^t')F(x,t')8(®(x,t')-v|/)G(t-0^' 

F(U 00 


Then we have 


J_ + ” p (x,t')U(x,t')f* (} ¥ ;x,0G(t-0dt , = ^VF U39 (y,w,x,t)dV = F 9 -{U(x,t)\vf) 

And the “complete” filtering should be 

C +oo / , i \ \ C "I" 1 "® C “I" 1 "® 

J (+o -(F(x,0|v(/))jv|/ = J ^ J ^ FF C/ 0 (F,V|/; x,t)dVd\\i 

f +00 s 

= J p(x,t f )U(x,t')G(t-t')dt' 

= p(j \c,t)U(x,t) 

Equation (46) can be extended to any other turbulent quantities, for example, VO, *S)( O): 


J_ + ” p(x,t')VO(x,0/i (r,x,f)G(t-t')dt' = F 9 -(vo|m/) 

P (*,*')$, (®(*,0)/i (w,x,t')G(t-t')dt' = F^ .(5 < (®)|v) = F <l> -5 f ( V ) 


(46) 


(47) 


(48) 


Where VO is viewed as a new random variable in addition to O. 

3.2 Transport Equation for Scalar DWFDF F 0 (v|/ ; x, t ) 

We can now derive the transport equation for scalar DWFDF from Equations (18), (19) as follows: 
first, we write the terms on the left hand side of Equation (19) as 


dp O m 

dt 


dr r dF 

-\y m F (S> dy=\y m -^-dy 


(49) 
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( 50 ) 


gp WK 

dx t 




II' I 1 ” 

—00 —00 

or 


V:—Fr 


SX ; 


u, o 


<7K di|/ 


CO / > 

K, 




t/V|/ 


Then, the terms on the right hand side of (19) can be written as 

' a “ 


_d_ 

SX: 


f <m) apO n A 


v 


8X: 


l J 


d 

dX : 




V 


= v„ 


J -GO 


8X ; 


i V 


p(m) 

dx, 


\d\\f 


i J 


(51) 


P S m = J S m (\v)F <]) d\\> = -^ 




g(^) 


<ii|/ 


(52) 


Where in Equation (52) we have applied the integration by parts and the following zero integration: 


—on ' 


if a finite value of filtering (A} exists. 


(53) 


Collecting all the integrand terms, and let the sum be zero (we refer this as the conservation condition, 
which is in general a sufficient but may not be necessary condition, Ref. 17), we obtain 


ggi, , g(Mvk)) _ f g U)gV 

dt dx { [fax dx i , 



k = l,2,---,M + 1 


(54) 


This equation also includes the equation of internal energy when S M+l (\|/) = 0 and other source terms in 
Equation (18) are neglected. Unlike the case of the joint DWFDF F U; 0 , for the marginal F® , the convection 
term is not closed because of the conditional filtering (u i | v|/^ . Then, this critically important term, 

corresponding to p U i Q> m in Equation (19), must be carefully modeled while the less important molecular 
diffusion term naturally remains in the closed form. In addition, we noticed that the equally important 
chemistry source term p S m (O) in Equation (19), which involves complex processes of turbulence- 
chemistry interaction, is closed in the scalar DWFDF equation, i.e., no need of modeling. This direct 
calculation of turbulence-chemistry interaction is one of the unique features of the PDF methodology. 


3.3 Nonlinear Model for Subscale Scalar Fluxes 

The convection term in Equation (54) contains the term • (Uj\\ |/) , which must be modeled. We may 

start from a nonlinear model, Equation (29), for the term p U i Q> m : 


p Ui® m =P Ui® m -T ( T m) 


0pO n 

SX: 


-T^-[c x Sij+c 2 nij) 


<3pO„ 

8X: 


This will lead to the following model by directly applying Equations (37) and (46): 


(55) 
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3.4 Summary 

With the model given by Equation (56), the scalar DWFDF equation for F 0 (i|/; x,t) can be written as 

(57) 

ifc = l,2,-»,Af + 1 

It can be verified that the scalar DWFDF Equation (57) can exactly deduce the time filtered equation (19). 
However, the model described by this equation is by no means unique. In addition, the variables ( U i , Sy 

and Clij ) are considered to be available during the solution procedure of the scalar DWFDF equation. 
Furthermore, in order to apply the available stochastic solution procedure built in the NCC code, we 
further simplify the model term in Equation (57) as follows: 





T 


where 


Fo dx i 


/V / ~ ^ 

l S ij ~\~ C 2 £2// 

p V 


c)x 

ux j ) 


SijS ij + Clij Clij 


Equation (59) is a crude approximation based on the dimensional argument for a time scale, which is 
responsible for the diffusion of scalar DWFDF in the sample variable space \| j k . We chose this time scale 
to be related to the rate of strain and rotation instead of its complex formulation. In order to prevent this 
time scale from being non-physically small during the simulation, we require that 

x > J(v + v T )/ £ (60 


because the right hand side of (60) represents the smallest time scale of the simulated flow field. 


4.0 Numerical Simulations of Single Element LDI Injector 

The lean direct injection (LDI) injector is a liquid fuel injector developed to reduce aircraft emissions. 
Stable combustion is essentially completed within a short distance through rapid fuel and air mixing. This 
design also allows for many small fuel injectors integrated into modules facilitating different fuel staging 
strategies, such as the one shown in Figure 1. So far, experimental observations have not fully clarified 
the dynamics of the mixing and combustion processes occurring in these injectors, and numerical studies 
need to be conducted to achieve a better understanding of the underlying physics of the LDI injector. 
Figure 2 shows the air s wirier and the convergent-divergent nozzle of the single element injector. Figure 3 
depicts the single element LDI combustor geometry and its computational domain. Five probes are 
dispatched along the centerline (Figure 4) to monitor the evolution of turbulent variables during the 
simulations. The numerical grid is formed using hexahedral elements and the total number of elements is 
about 862,000, which is a relatively coarse grid used in a previous RANS simulation (Ref. 18). 
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In this study, the liquid fuel is Jet-A, and C12H23 is adopted as its surrogate, the fuel is injected at the 
throat of the nozzle, mixing with the swirling air that comes from the air swirlers which consists of six 
helical, axial vanes with downstream vane angles of 60°. A prescribed droplet-size distribution spray 
model is used. 

In this section, we present the VLES results of TFNS approach with and without invoking Eulerian 
scalar DWFDF method. We started both simulations with a same initial flow field that was created by 
earlier simulations from URANS to TFNS at the end of time step of 93,000. In the pure TFNS 
simulations, the chemistry-turbulence interaction is modeled with a “laminar chemistry” formulation, in 
which the generation rate of compositions is determined by the “known’ filtered turbulent variables. This 
simulation has been carried out up to 185,000 time steps with a time step dt = lxlO 6 . In the TFNS 
simulations with invoking an Eulerian scalar DWFDF method, the compositions and internal energy will 
be determined by the scalar DWFDF equation while the flow field is still determined by the continuity 
and momentum equations of the filtered Navier-Stokes equations. Since the production rate of 
composition (chemistry-turbulence interaction) is in a closed form in the scalar DWFDF equation, it can 
be directly calculated without modeling. This simulation is carried out up to 198,000 time steps with the 
same size of time step. The stochastic numerical procedures for solving Eulerian scalar DWFDF equation 
and fuel spray equation are described in Reference 4 and adopted in the present simulations. 

Figure 5 shows the general pictures of simulated spray reacting flow in a single element LDI injector 
by both methods. Most of the results will be presented side by side for comparisons. The main results are 
presented in terms of : 1) the convergence history of time accurate simulations, which reflects numerical 



Figure 1. — Nine element LDI module. 


Figure 2. — Single element LDI injector. 




Figure 3. — Computational domain and grid spacing. 


Probe locations: 

1 . x = 0.00450m or -2.7 mm from dump face 

2. x = 0.01015m or 3 mm from dump face 

3. x = 0.0191 5m or 1 2 mm from dump face 

4. x = 0.13150m or 124 mm from dump face 

5. x = 0.18715m or 180 mm from dump face 


12 3 45 



Figure 4. — Probes located along centerline. 
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Figure 5. — Global pictures of spray reacting flow simulated with and without scalar DWFDF method. 
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Figure 6. — Convergence histories of time accurate simulations. 


performance, 2) the variations of velocity components and temperature versus the time step at five probes 
along the centerline (see Figure 4), which indicate the development of the simulated reacting flow, 3) the 
centerline distributions of mean (i.e., time averaged) temperature and mean axial velocity, 4) the mean 
temperature and mean velocity profiles along Z (crossing) axis at several downstream locations, 5) the 
center recirculation zone visualized by an iso-surface of zero mean axial velocity and the contour plots of 
various mean and instantaneous turbulent quantities in the X-Z center plane. These profiles, iso-surface 
and contour plots will provide the information about flow and flame structures of the simulated spray 
reacting flow. Some available experimental data (Ref. 19) are also plotted for comparisons with the 
numerical simulations. 


4.1 Convergence History 

Figure 6 shows the time history of the number of subiterations incurred in time accurate simulations 
with and without invoking the Eulerian scalar DWFDF method over the last 15,000 time steps. It is 
interesting to see that the simulation with the Eulerian scalar DWFDF method converges faster than the 
simulation without invoking scalar DWFDF method. For each time step, the former converges after about 
20 subiterations, but the later often needs much more subiterations to converge. Comparing the wall time 
of computing, the simulation with DWFDF is about 34 hr with 256 processors for the 15,000 time steps; 
however, the simulation without DWFDF needs over 40 hr to finish 15,000 time steps of calculations. 
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4.2 Time Variation of Instantaneous Variables at Centerline Probes 

Time history of temperature and velocity components at five downstream locations are recorded during 
the simulations (see Figure 7), it can be seen that at the first three locations the flow appears to be strongly 
fluctuated and fully developed turbulence (see both the instantaneous temperature T and axial velocity u ). 
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Figure 7. — Time history of temperature T and velocity components u, v and w at 5 probes. 
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Figure 7. — Concluded. 
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4.3 Centerline Distributions of Mean Axial Velocity and Temperature 

Centerline distributions of mean axial velocity and temperature are shown in Figure 8 and Figure 9. 
The effect of scalar DWFDF equation on the simulation is noticeable and positive. 




Figure 8. — Centerline distribution of mean axial velocity U. 
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Figure 9. — Centerline distribution of mean temperature T. 
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4.4 Mean Temperature and Velocity Profiles Along Z Axis at Downstream Locations 

Mean temperature profiles at downstream locations at x = 5, 10, 20, 50, 110 and 200 mm are 
compared and shown in Figure 10. These comparisons have revealed a positive effect of the DWFDF 
equation on the prediction of temperature. However, the effect of DWFDF equation on the velocity field, 
shown in Figure 1 1 and Figure 12, is mixed. 



Figure 10. — Comparison of temperature distribution along Z axis at downstream locations. 
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Figure 10. — Concluded. 
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Figure 1 1 . — Comparison of velocity component V along Z axis at downstream locations. 
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Figure 11. — Continued. 
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Figure 11. — Continued. 
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Figure 11. — Concluded. 





Figure 12. — Comparison of velocity component U along Z axis at downstream locations. 
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Figure 12. — Continued. 
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Figure 12. — Continued. 
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4.5 Flow Structure and Contour Plots of Variables in the Center X-Z Plane 

To reveal the simulated mean flow and flame structures, we have plotted the center recirculation zone 
using the iso-surface of zero axial velocity (see Figure 13), the contour plots of velocity components, 
vorticity magnitude, temperature, turbulent kinetic energy and species mass fractions Ci 2 H 2 3 , 0 2 , C0 2 , 
which are shown in Figure 14, Figure 15, and Figure 16. From these figures we observe that the structures 
of all scalar fields simulated by TFNS with invoking DWFDF equation are quite different from the ones 
without invoking DWFDF equation. Also, we notice that the level of subscale turbulent kinetic energy is 
significantly reduced in the TFNS simulation with invoking DWFDF equation (a phenomenon often 
found in the simulations with the conventional or standard FDF method); however, it does not severely 
affect the global structure of turbulent reacting flow in our simulations due to the adoption of the new 
DWFDF equation, in which the time scale or frequency for the DWFDF diffusion in the sample space is 
not directly related to the subscale turbulent kinetic energy and its dissipation rate, instead, it is 
determined by the strain and rotation rate of the large turbulent flow, see Equation (59). The contours of 
all instantaneous variables in X-Z plane are plotted in Figure 17. 
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Figure 13. — Comparison of center recirculation zone. 
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Figure 14. — Comparison of mean velocities U, V, W and vorticity contours in center plane. 
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Figure 14. — Concluded 
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Figure 15. — Comparison of mean temperature and k contours in center plane. 
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Figure 16 . — Comparison of C12H23, O2 and CO2 contours in center plane. 
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Figure 17. — Contours of instantaneous variables in X-Z center plane. 
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Figure 17. — Continued. 
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Figure 17. — Concluded. 


5.0 Conclusion 

The very large eddy simulations of Jet-A spray reacting flow in a single element LDI injector have 
been carried out both with and without invoking an Eulerian scalar DWFDF method. The NCC code is 
used under the same numerical setting, same computational geometry and grid resolution. In addition, the 
same spray model and the same chemistry kinetics are used for all the simulations. In this way, we hope 
to isolate and study only the effect of Eulerian scalar DWFDF method on the simulations, paying a 
particular attention to the simulation quality and numerical performance. 

From the present simulations we have observed that there are noticeable improvements on the 
temperature prediction in the region of strong turbulence by using the Eulerian scalar DWFDF method. 
The distribution or structure of other scalar quantities appears to have a significant difference between the 
simulations with and without the scalar DWFDF method, and it seems to have more reasonable results by 
invoking the scalar DWFDF method. However, we need more reliable experimental data to evaluate and 
draw the conclusion. For the flow field, the effect of the scalar DWFDF method is mixed. Better 
predictions were observed in the axial velocity component. 

As to the numerical performance, the observation is that the simulations with invoking an Eulerian 
scalar DWFDF method appear to have a faster convergence and more stable than the simulations without 
invoking Eulerian scalar DWFDF method. 
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